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ABSTRACT 

The temperature of the intergalactic medium (IGM) is set by the competition between photoheating 
and adiabatic cooUng, which are usually assumed to define a tight equation of state in which the 
temperature increases monotonically with density. We use a semi-analytic model, accurate at low 
to moderate IGM densities (< 5 times the mean), to show that inhomogeneous reionization can 
substantially modify these expectations. Because reionization is driven by biased sources, dense 
pockets of the IGM are likely to be ionized first. As a result, voids initially remain cool while dense 
regions are heated substantially. However, near the end of reionization, dense regions have already 
cooled from their initially large temperature while voids have only just been heated. Thus, near the 
end of helium reionization the equation of state can invert itself, with the hottest gas inside voids. 
The degree to which this happens depends on the magnitude of the density bias during reionization: if 
rare, bright sources dominate reionization, so that each ionized region contains a typical volume of the 
IGM, the equation of state will remain monotonic. We also show that the distribution of temperatures 
at a fixed density has significant scatter and evolves rapidly throughout and even after reionization. 
Finally, we show that the observed temperature jump at z ~ 3.2 is consistent with the behavior at 
the end of helium reionization, although it requires a somewhat larger temperature increase than 
expected. 

Subject headings: cosmology: theory - intergalactic medium - diffuse radiation 



1. INTRODUCTION 

The two most dramatic events in the history of the 
intergalactic medium (IGM) are the reionization of hy- 
drogen (by the first generations of galaxies) and he- 
lium (by quasars). These have become key landmarks 
for both observational and theoretical cosmologists in 
the past several years. Evidence for hydrogen reioniza- 
tion comes from a number of directions, none of them 
clear but all consistent with (p ossibly extended) reion- 
ization at z ~ 6-10 fsee lFan et al. 2006; Furlanctt o et al.l 
l2006bl for recent reviews). The clearest evidence for 
helium reionization comes from He II Lya forest spec- 
tra in the extreme ultraviolet, which seem to show a 
rapid evolution in the characteristics of the transmis- 
sion at z ~ 3 JPavidscn ct al. 1996; Anderson et al. 199^; 



Heap et al . 2000; S mette et a l. 2002; Zheng et al. 2004b; 



Shull et al.ii2004:.Reimers et al...2004, : .Zheng et al..,2004a : 
Reimers et al.ll2005D . although there is also some indirect 



evidence (see below). 

These two reionization epochs are largely responsi- 
ble for determining the thermal history of the IGM. 
Before hydrogen reionization, the neutral IGM cooled 
adiabatically until the first structures formed, proba- 
bly reaching temperatures T ~ 30 K. X-rays from 
the first galaxies most likely slowly he ated the neutral 
IGM, to T < lOOOK (p ^h 2001; VcnkatcsanetaDIlOOl 
iKuhlen et all 120061 : iFurlanctto 2006). However, hydro- 
gen reionization was a much more dramatic change: the 
^ 1 eV leftover from ea ch ionizing photon heated the 
IGM to ~ 1-2 X lO'^ K (|Miralda-Escude fc ReesI [l99l 
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I Abel fc Haehlleltl [l99l pTittlev fc Meiksh] [2006h . The 
harder photons responsible for helium reionization would 
have rehe ated the IGM to sim ilar, or even larger, tem- 
peratures (|Hui fc Gnedinlll997i) . 

Once reionization is complete, this heating channel 
slows dramatically - because only the relatively small 
fraction of ions that recombine will couple to the pho- 
toionizing background. The subsequent temperature 
evolution is determined by a combination of adiabatic 
heating and/or cooling, photo-h e ating , and Compton 
coolin g (|Miralda-Escude fc Reei Il994t iHui fc GnedinI 
119971 ).'^ The competition between these processes forces 
the gas temperature to approach an asym ptotic form 
set by the background ion izing spectrum (jHui fc GnedinI 
fTMTtlHuTfc Haimanl[2003l) . Because the magnitude (and 
indeed sign) of the adiabatic term depends on whether 
the gas is over- or underdense, the IGM assumes an equa- 
tion of state T K, Tq{1 + (5)'''"^, where 5 is the fractional 
overdensity of the gas element and 7 is nearly indepen- 
dent of 5 in most models. 

This equation of state inevitably affects many observ- 
ables of the Lya forest, including th e power spectrum 
of the transmitted flu x (Croft 1998; Zaldarr iaga et all 
[2OOII: iViel et all [200l [McDonald et al. 20061 . the op- 
tical depth distribution ( Lidz et al.l 120061), and the line 
widths of i ndividua l features (ISc have et al.l Il999l 120001 : 
iRicotti et a l. 2000; McDonald et al. 2001). In particular, 
the observed temperature evolution has been used to con- 
strain the epochs of helium and hydrogen reionization: 
most dramatically, the velocity widths of L ya forest ab- 
sorbers seem to increase sharply at z '--^ 3.2 ( Schave et al.l 



More precisely, several other factors contribute (see the Ap- 
pendix of Mui & Gncdin 1993 for ^ summary) but these three are 
by far the most important. 
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l200Cir). while the equation of state simul taneously flattens 
(|Schave et ani2000l : [Ricotti et al.ll2000l ). This may be a 
result of helium reionization. If so, the temperature at 
higher redshifts could provide a clean "fossil" record of 
hydrogen reionization: it appears to be relatively large 
and so requires that epoch to have occur ed at 2; < 10 
(|Theuns et al.ll2002cl : [Hui fc Haimanll2003f) . 

The cartoon history described above has been devel- 
oped from simple uniform reionization models. These ex- 
isting treatments have ignored one extremely important 
facet of reionization: it is highly inhomogeneous, espe- 
cially bec ause the sources driving it are themselves highly 
clustered (i Barkana fc L oebl 120041 : iFurlanetto et"ani2004l : 
iFurlanetto fc Oh 2007). As a result, overdense regions 
in the IGM are ionized before the empty voids (at least 
on average), and their thermal histories will evolve dif- 
ferently. This affects the mean equation of state and also 
adds an element of "stochasticity" to it, because patches 
at identical densities can have different reionization histo- 
ries and hence different temperatures. To this point, an- 
alytic models have ignored this aspect for simplicity, and 
numerical simulations have not been able to include both 
the inhomogeneity of reionization (which occurs on large 
> 10 Mpc scales) while also resolving features of the Lya 
forest. The exception is the M onte Carlo, semi-analytic 
helium reionization model of Gles er et al.l |2005), who 
showed that this inhomogeneity significantly affects the 
equation of state. However, they used a highly approxi- 
mate prescription for inhomogeneous reionization. 

In this paper, we present a model that incorporates the 
latest models of inhomogeneous reionization for calculat- 
ing the IGM equation of state. We apply this model to 
the helium reionization epoch and show that it has an 
enormous impact on the thermal properties of the IGM 
during the best-observed epochs, z ~ 2-4. We describe 
our model for inhomogeneous reionization in ^}2] and our 
model for the subsequent temperature evolution in f}3l 
We examine the resulting temperature distributions in 
21 the equation of state in fj5l and the evolution of the 
IGM temperature in S}6l We compare our models to ob- 
servations in f}7] and conclude in fj8l 

In our numerical calculations, we assume a cosmol- 
ogy with = 0.26, = 0.74, = 0.044, H = 
lOO/i km s"^ Mpc-i (with h = 0.74), n = 0.95, and 
as = 0.8, consistent with the most recent measurements 
(jSpergel et al.ll2007t ). Unless otherwise specified, we use 
comoving units for all distances. 

2. INHOMOGENEOUS REIONIZATION 

Our first goal is to compute the distribution of reion- 
ization redshifts for gas elements with mass rrip at a 
linearized fractional overdensity, extrapolated to the 
present day, of S^. For notational simplicity, we will use 
S = a'^{m), the variance of the density field smoothed on 
mass scale m, as a proxy for mass in this section; also, 
Sp = S{mp). We will include both helium and hydrogen 
reionization; for clarity, we will refer to the moment at 
which ionized hydrogen fills the universe as zu and the 
moment at which doubly-ionized helium fills the universe 
as zhc-'' 

For most of our calcula tions, we will use the ex cursion 
set reionization model of IFurlanetto et al] (|2004[ ). which 

* We assume that hcUum is singly-ionized at zjj as well. 



provides an excellent match to simulations of hydrogen 
reionization, at least when IGM recombinations are rela- 
tively uniform (|Zahn et al.ll2007l: IMesinger fc Furlanettcl 
l2007f ). According to this model, under the reasonable 
assumption that ionizing photons are produced inside 
galaxies, the ionization field can be built from the lin- 
earized density field by demanding that a region is ion- 
ized if and only if it has 



CfcolliSl^S) >l 



(1) 



where ^ is the number of ionizing photons escaping into 
the IGM per atom (of the appropriate element) inside 
of galaxies,^ fcon{Sj^,S) is the collapse fraction in a re- 
gion of density 61 and mass S, and N^ec is the aver- 
age cumulative number of recombinations per ionized 
atorn . Using t he excursion set formalism (jBond et al.l 
II99II : iLacev fc~ Colc 1993|), this condition can be trans- 
formed into the comoving number density of ionized bub- 
bles with masses in the range (to, to^ -|- dm): 



m{m,z) = \ll ^ 



dine 



dlni 



Bojz) 
a{m) 



exp 



B^{m, z) 

' 20-2 (to) 



(2) 

where the excursion set barrier has been approximated 
by the linear function B{m,z) « Bf){z) + Bi{z)S{m) and 
p is the mean comoving mass density. Note that, unlike 
the spherical collapse barrier, this B{z) only extends to 
a minimum mass scale TOb^min equal to C, times the mini- 
mum galaxy mass mmin; smaller ionized bubbles cannot 
contain any ionizing sources so do not appear in the for- 
malism. 

This reionization model has been applied extensively to 
the epoch of hydrogen reionization, and we will also use 
it for helium reionization. There are some key differences 
that make the model less appropr iate for the helium era, 
but it is nevertheless useful (see IFurlanetto fcOhll2007l 
for a closer look at these issues). First, the excursion set 
formalism divides the universe into completely ionized 
and completely neutral phases. This is an excellent ap- 
proximation for hydrogen reionization, where the mean 
free paths of ionizing photons are much smaller than the 
ionized bubbles (assuming that ultraviolet photons drive 
the process). It is less good for helium reionization: be- 
cause quasars have hard spectra, many of the ionizing 
photons are able to travel large distances (many Mpc) 
before interacting wtih a Hell ion. Thus the boundaries 
of ionized bubbles will be somewhat "fuzzy" during he- 
lium reionization. However, the mean free path of most 
ionizing photons is considerably smaller than the scale 
of the ionized features throughout most of reionization 
(e.g., a typical quasar can ionize a region ~ 10 Mpc 
across, while ~ 50% of ionizing photons are absorbed 
within ^ 3 Mpc). 

Another difference is that helium reionization is mainly 
driven by relatively rare, bright sources. As a result, ran- 
dom fluctuations in their abundances play an important 
role, and because their helium bubbles are so large they 
will contain a more random mix of dense gas and under- 
dense voids than the bubbles driven by small, clustered 
sources during hydrogen reionization. Both these factors 

^ For reference, our "slow" and "fast" helium reionization models 
require fjjo = 10.2 and 33.2 if zhc = 3, and hydrogen reionization 
requires C^h = 40 if = 8. 
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imply more stochasticity in the reionization history for 
a given patch density than in our density-driven model. 
In practice, the early stages of helium reionization are 
primarily stochastic, but once the ionized fraction ex- 
ceeds one-half or so, the density-driven component takes 
over. We will therefore contrast it with a model in which 
reionization is completely independent of the local den- 
sity (but in which the mean ionized fraction evolves the 
same way). 

This distribution allows us to compute the range of ion- 
ization histories for our pixels through Bayes' theorem, 
which in this case says that 

IZ ■ '^rnirn/p)nb{m,z)ps^\s{Sl\B) 



PS, 



(3) 

Here Ps^ {zicm > z\d1) is the probability that the element, 
with at mass Sp, was ionized before z, pSpi^l) is the 
probability that an element of mass Sp has density 
(this is simply a gaussian with variance Sp), and 




Fig. 1. — Cumulative and differential probability distributions 
of the helium reionization redshift. We assume that helium reion- 
ization completes at zjje = 3. In both panels, the thick curves 



^2tt{Sp - S) 



correspond to mass elements with 59 



-2.5, 0, 2.5, and 5, 



X exp 



(SI 



B) 



2 1 



2{Sp - S) 



from left to right. In the top panel, the thin curve shows the 
globally-averaged reionization history for our fiducial model. 



(4) 



is the conditional probability that an element with mass 
Sp has density (5°, given that it has density B at mass 
S < Sp. In other words, the cumulative probability that 
an element has been ionized before z is the sum of the 
probabilities that it lies on trajectories passing through 
the excursion set barrier B(z), normalized by the total 
abundance of elements of the appropriate final density. 
Equation fS]) can be simplified as 



PsMon>z\Sl)^ 



Bn 



dS 



St) 



V S^ 



■ exp 



{5S 



Here Sb = <t (wb^min) is the variance of the density field 
at the mass scale corresponding to the smallest allowed 
ionized bubble. Obviously, if Sb > Sp we must have all 
the pixels with 51 > B ionized. In reality these scales 
are not too far apart: mb^rnin is C times the minimum 
galaxy mass at the appropriate redshift, while is at 
the Jeans scale for ionized gas. Their ratio is: 

1 o I ^ \ f ^vir.mm \^^^ ( TiGM , 

1-3 I ) ( o inRT^ ) ( ) ' (6) 



2SpS{Sp 
/-N tions 
(5) 



105 K 



10* K 



where T,, 



is the minimum virial temperature of ion- 



izing sources (the fiducial value corresponds to the mini- 
mum galaxy mass in photoheated gas). The similarity 
between these two scales should come as no surprise, 
since mb^min is definined as the minimum mass which 
can accrete in a photo-ionized IGM, which is of course 
intimately related to the Jeans mass. This leads to some 
ambiguity in our results (see below). 

We show the resulting distributions of the local helium 
reionization redshift for several different in Figure [T] 
(this era has much more significant implications for ob- 
servations than hydrogen reionization); the upper and 
lower panels show the cumulative and differential proba- 
bility distributions, respectively. We have assumed that 



reionization ends at zjjc = 3. In each panel, the thick 
solid, long-dashed, short-dashed, dot-dashed, and dotted 
curves are for elements with = —5, —2.5, 0, 2.5, and 
5, respectively. We take nip to be the Jeans mass of gas 
in an ionized medium of temperature 10^ K at mean den- 
sity for z = 3; note that a{mp) w 2.85 at this epoch, so 
these are approximately one- and two-sigma fluctuations 
at the redshift of helium reionization. 
The inhomogeneity of reionization is obvious from Fig- 
BS^^ ] even at a fixed density, gas elements can have a 
r'ange of ionization histories. However, the distribu- 
are strongly dependent on the underlying density, 
with denser gas much more likely to be ionized at high 
redshift. This is simply because such regions are usually 
found in large-scale overdensities, which are ionized first 
so long as ionizing sources are more biased than recom- 
bining clumps. 

2.1. Dependence on the Reionization Model 

The thin curve in the top panel of Figure [T] shows the 
globally-averaged reionization history, Xi{z). It mirrors 
the distribution of reionization redshifts of mean den- 
sity pockets for much of the evolution, although it has 
a significantly longer tail (to account for dense gas that 
is ionized early on). If zhc is fixed, the duration of the 
reionization era is determined by the rate of change of 
/coll. In our fiducial model, we allow all dark matter 
halos with m > mmin = 'm-i to form quasars, where rrii 
is the minimum mass for dark matter halos to accrete 
material in a photoionized medium (with T ^ 10"* K); it 
corresponds to Tvir.min ~ 2 x 10^ K. 

This is roughly consistent with other h elium reioniza- 
tion models (e.g., IWvithe fc Loebl l2003l ). but reioniza- 
tion may be less extended if more massive halos domi- 
nate the ionizing photon budget. This can occur in two 
ways: by increasing the efficiency of more massive galax- 
ies or by increasing m,^i^. As an example, in Figure [2] we 
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Fig. 2. — As Fig. [T] except rrimin = lOrdi. 



show a model with mmin = lOm^. In this case, ~ 70% 
of ionizing photons are produced after = 4 (assum- 
ing zhc = 3). This provides a somewhat better match 
to the reionization history inferred from observations of 
quasar abu ndances (in particula r, the quasar luminosity 
function of iHopkins et al.ll2007f ). although uncertainties 
in the faint end are rather large at high redshifts. Al- 
though the reionization history is faster, note that the 
qualitative features are the same, and the dependence 
on density remains substantial. 

As described above, another key component of a he- 
lium reionization model is the degree of stochasticity in 
the process. The models described above do not include 
this component, but it is easy to construct such a model: 
for every density, we simply draw P{> z) from the distri- 
bution of Xi{z). Thus, the histories of every gas element 
in these models tend to mimic those at the mean den- 
sity of the biased models, but with longer tails toward 
higher redshift. We will compare these kinds of models 
extensively below. 

2.2. Dependence on the IGM Jeans Mass 

Figure [3] illustrates the ambiguity in choosing rup. We 
show three sets of curves, with 5\ = —5, 0, and 5, from 
left to right. Within each set, the solid, dashed, and dot- 
ted curves take T = 2 X 10^, 10^, and 10^ K when evalu- 
ating the Jeans mass of the IGM.^ The temperature has 
the largest effect on high-density regions: smaller Jeans 
masses give more "space" between Sb and S'p, allowing 
trajectories to wander far from the ionization barrier be- 
yond its low-mass edge. This allows high-density regions 
to be ionized later than otherwise (by trajectories that 
begin at low densities and then wander to high densities 
only at small masses). 

® These three choices bracket the expected post-reionization 
IGM temperature (as illustrated in ^ below): a high temper- 
ature would be possible with recent hydrogen reionization and 
a substantial X-ray background; the 10'^ K choice is still below 
the temperature exp ected with very early hydrogen reionization 
IIHui fc Haimanll200a') . 
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Fig. 3. — Cumulative probability distributions of the helium 
reionization redshift. We assume that helium reionization com- 
pletes at = 3. The three sets of curves correspond to mass 
elements with 5^ = —5, 0, and 5, from left to right. Within each 
set, the solid, dashed, and dotted curves take T = 2 X 10"', 10"', 
and 10^ K when evaluating the Jeans mass of the IGM. 



The trend is reversed for low-density pixels: here a 
smaller Jeans mass translates to somewhat earlier reion- 
ization. Again, this can be understood in terms of the 
diffusion problem: if Sp 3> Sb^ trajectories have an en- 
hanced opportunity to wander from high densities (where 
they crossed the barrier) to the low-density pixel. How- 
ever, the overall difference is smaller in voids than in 
dense regions, because trajectories must still wander 
quite far from the ionization barrier to end up as voids 
at the pixel scale (whereas they can begin just below the 
ionization barrier and easily wander to high densities). 
The two effects roughly cancel at average densities. For- 
tunately, this ambiguity does not have a major impact on 
our final results, which are much more sensitive to the 
thermal evolution of voids and average-density regions 
than to overdensities. 

3. THE TEMPERATURE HISTORY 

Now that we have obtained the distribution of ioniza- 
tion histories, we wish to compute the subsequent tem- 
perature evolu t ion. W e will follow a simplified version of 
IHui fc Gnedil] (fl99l . 

Consider a gas element of fractional overdensity 5 (note 
that this is the true nonlinear overdensity) and temper- 
ature T illuminated by an ionizing background of the 
form 

J. = JHI.-21 f — X ( 1 < ^Hell (7) 

where Jhi,-21 is the angle-averaged specific intensity of 
the background, in units of 10~^^ erg cm~^ s~^ Hz"^ 
sr^^, uai and t'Heii are the frequencies corresponding 
to the HI and Hell ionization edges, and / is a con- 
stant. We will assume for simplicity that / = be- 
fore helium is reionized and / = 1 afterwards. We 
therefore do not evolve the magnitude of the ionizing 
background with redshift; our results are independent of 
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the a mphtude of the back ground in the highly-ionized 
hmit (jHui fc GnedinI fT997h . although our simple treat- 
ment will not be completely accurate if the shape of the 
spectrum evolves substantially (which it probably will as 
quasars begin to dominate over galaxies). However, this 
only has a small effect on the calculated thermal history. 
Our assumed radiation background in equation ([7]) only 
affects the species abundance via the assumption of pho- 
toionization equilibrium; the heating rates - which are 
a complicated function of radiative transfer effects - are 
determined separately. 

For simplicity, we will take a fiducial spectrum with 
a = 1.5 and Jhi,-21 — 1, similar to that expected from 
quasars. This yields photoionization rates, defined via 



' 47rJ, 
ai' cr," 



(8) 



of rHi,-i2 = 2.82 (2.83), rHoi,-i2 = 4.03 (4.29), and 
rHcii,-i2 = 0(0.092) before (after) He II reionization, 
where F — Fi__i2 x 10~^^ s~^. This spectrum is slightly 
harder tha n recent estimates of quasars in the far- 
ultraviolet (jVanden Berk et al.|[200lt iTelfer et an[200l 
{a ~ 1.6-1.8, albeit with a wide dispersion), but we 
use it for consistency with past estimates; spectra in 
the above range do not appreciably affect our results, 
although they do slightly decrease the asymptotic tem- 
peratures. Note also that, before quasars dominate the 
ionizing background, this power law is much too hard for 
hydrogen ionizing photons, so we do not expect to get the 
early temperature history precisely correct. However, we 
focus on the era around helium reionization, when it is a 
good approximation. 

We define the number density of species i to be = 
(1 + S)XiPb/mp, where pt is the proper mass density of 
baryons. Note that the Xi are the species fractions by 
mass relative to all baryons (so Xm is not the neutral 
fraction). 

The thermal evolution of each element is determined 
by 



2T dS T d(E^l, 



^ = -2HT . 

dt 3(1 + S) dt Y.. X, 



dQ 



dt 



Sksnt dt ' 
(9) 

where d/dt is the Lagrangian derivative. On the right 
hand side, the first term accounts for the Hubble expan- 
sion, the second for adiabatic cooling or heating from 
structure formation, the third for the internal energy 
gained or lost per particle from changing the total par- 
ticle density, and in the last term dQ/dt is the net heat 
gain or loss per unit volume from radiation processes (see 
below) and nt is the local (total) particle number density. 

The second term requires a n expression for the g rowth 
of the nonlinear density field. iHui fc GnedinI l)1997l ) used 
the Zeldovich approximation, along with an analytic es- 
timate for the distribution of the strain tensor, which is 
probably the best analytic approximation to the density 
evolution in the quasilinear regime. However, to main- 
tain computational simplicity, and so as not to introduce 
another probabilistic clement to our code, we will assume 
that the gas elements evolve following the spherical col- 
lapse (or expansion) model. We map linear densities to 
nonlinear overdensities via the following fitting formula 



due to IMo fc Whiti (fl99l . 

. 1-35 1.12431 



0.78785 



■ (1 + (5)2/3 (1 + 5)0-58661' 

(10) 

where D{z) is the linear growth function and Sc ~ 1.69 
is the threshold for virialization in the spherical collapse 
model. This fitting formula fails at shell-crossing, but 
we will not be considering such extreme densities. Under 
this assumption, we can write 



dS _ . dD dS 
dt " ° 'dtdS2' 



(11) 



The fourth term in equation Q includes a number 
of radiative heating and cooling processes. The most 
important heating mechanism is photoionization itself; 
each species i contributes a term 



dQ 
dt 



di/ {A-Kj^)(Ji 



(12) 



where ai is the photoionization cross section for species 
i, Vi is its ionization threshold, and a Roman h denotes 
Planck's constant (to differ entiate it from the Hubble 
constant). We use the fits of IVerner et al.l ()1996[ ) for the 
photoionization cross sections. The other relevant mech- 
anisms cool the gas and include Compton cooling off the 
CMB, recombinations (radiative and, for He II, dielec- 
tronic), coUisional ionization, collisional line excitation, 
and fr ee- free emission. We use the fits of iHui fc GnedinI 
(|1997f ) for all of these except Compton cooling (for which 
we use the exact form) and free-free emission (for which 
we use the fit in iTheuns et al.l [1998*) . We have verified 
th at our results are unc hanged if we use the fit s presented 
in ITheuns et~an (|1998l updated from those in lCenlll992D 
for all the listed mechanisms. 

Because the last two terms in equation ([9|) depend on 
the species abundances, we must supplement them with 
equations for the evolution of each Xi . These are partic- 
ularly important during the early stages of reionization, 
when the Xi change rapidly, along with the photoheat- 
ing and radiative cooling rates. However, after this initial 
phase the abundances rapidly settle into photoionization 
equilibrium and subsequently evolve quasistatically with 
the Hubble expansion. Unfortunately, the thermal evo- 
lution during reionization depends on radiative transfer 
effects and is not well-described by o ur simple model. 
Thes e tend to harden the spectrum llAbel fc Haehnel'3 
I1999D . because gas near the ionizing sources preferen- 
tially absorbs soft photons, leaving more distant gas to 
be ionized by high-energy photons. Our procedure in- 
stead treats each element in isolation. 

For simplicity, rather than solve the dXi/dt equations 
explicitly, we initialize our models shortly after reioniza- 
tion by assuming photoionization equilibrium at a fixed 
temperature T^. We then assume that the gas remains 
in photoionization equilibrium at each temperature, den- 
sity, and redshift. Beyond the initial phase, this provides 
excellent agreement with full calculations including the 
species rate equations but is considerably faster to com- 
pute. 

We must then specify Ti for both HI and Hell. For 
concreteness, we take = 2.4 x 10* K throughout; this 
is a rather large value (implicitly assuming significant 
hardening from radiative transfer, or intrinsically hard 
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reion ization sources such as quasars; lAbel fc Haeh^ 
I1999D . However, our results are quite insensitive to it, 
because the subsequent cooling is so rapid that memory 
of the initial post-reionization temperature is erased well 
before the z ~ 2-4 regime relevant for us (T heuns et al.l 
l2002cl : lHui fcHaiman 200a ). Indeed, we ran our models 
for smaller values of and found little difference in the 
results, except that the pre-helium reionization temper- 
ature was slightly smaller (and therefore disagreed even 
more with measurements; see Sj7]). 

The choice of T^'^ is much more important, because 
it directly impacts the temperature in the observable 
redshift range. Some general arguments do provide 
bounds on its likely value (jMiralda-Escudc & Roes 1994; 
lAbel fc Haehnelti 11999). If the region between the gas 
element and the ionizing sources is optically thin, the 
mean excess energy of ionizing photons £'thin is the av- 
erage of the entire spectrum weighted by Oi oc ^ so 
(-Ethin) lEi K, {a + 2)~^ « 0.3, where Ei is the ioniza- 
tion potential. If instead we take the optically thick 
limit, where all ionizing photons are absorbed, the pho- 
ton energies are unweighted such that (-Ethick) lEi k, 
(oL — V)~^ ~ 2. In either case, this energy must then 
be shared with all the IGM baryons through Coulomb 
interactions (recall that only « 0.07 of electrons begin 
inside of Hell atoms); the net temperature change would 
then be AT w 0.035(2/3/cb) {E) - 0.5-3 x lO"* K with 
our assumed spectrum. 

We will choose = 3 x 10^ K for our fiducial model 
(corresponding to AT '--^ 2 x 10"* K after accounting for 
the initial temperature), but we will also explore the ef- 
fects of larger and smaller values. Note also that we 
treat T^'^ as independent of the underlying density, even 
though the initial temperature is a (weak) function of the 
underlying density and AT is actually the constant for a 
specified spectrum. We have made this tradeoff in the in- 
terests of computational simplicity and transparency. In 
fact, we expect both the amplitude and shape of the ion- 
izing background to be density dependent; the radiation 
field both weakens and hardens as it pr opagates outwards 
from the source (e.g.. iBolton et al][2004) . To the extent 
that voids lie farther from the ionizing sources than dense 
filaments, they will receive more highly-filtered radiation 
and hence have larger AT. This hardening effect, which 
we ignore, amplifies the effect we now describe: an in- 
verted equation of state in which voids are substantially 
hotter than in a homogeneous reionization scenario. On 
the other hand, at z < 3 there is some evidence for a 
softer radiation field in voids than i n filaments, whic h 
would weaken the effect we describe (Sh ull et al.ll2004|) . 

We must note that equation ([9]) neglects some phys- 
ical mechanisms that may affect the thermal evolution 
of certain gas particles. Most important is shock heat- 
ing, which dominates the thermal budget of the IGM at 
low r edshifts [z < 1; ICen fc Ostrikei1ll999l : iDave et all 
1200 It ) and whic h cannot be enti r ely ignored even at 
higher redshifts (iNath fc Silkl l2001l: i Valageas et a"n i2002t 
iFurlanetto fc Loebl 12004^ 7" Detailed simulations of large- 
scale structure shocks show that they typically sur- 
round sheets, filaments, and halos but can also extend 
into low density gas on rare occasions (|Rvu et all 120031 : 
iKang et al.ll2005l ). 

However, such shocks are likely to have only a small 



effect on the Lya forest (and diffuse IGM) at z > 2, 
because nearly all of the forest - with the exception of 
rare, high column density systems - has densities near or 
even below the mean. In this regime, even analytic mod- 
els that assume photoionization equilibrium and explic- 
itly ignore sho cks have had a good deal of success (e.g., 
IBi et al.lflMl lBi| [T99lls"chavJ200iri The reason can be 
seen in, e.g.. Figure 11 of lDave et al. (|1999D . which shows 
the IGM equation of state measured in a cosmological 
simulation at a variety of redshifts; this includes large- 
scale structure shocks but excludes helium reionization, 
of course. At z = 3, the vast majority of gas is confined 
to the well-defined equation of state. Shocks, which raise 
the temperature above that expected from this relation, 
affect only a small fraction of gas, and virtually all of it 
has A>10. Atz = 2, more of the gas is in the shocked 
phase, but again virtually all has A > 10. That shocks 
only affect gas sitting at higher densities than we will 
study here is not entirely coincidental: naively, shocks 
should occur when flows begin converging. In the cos- 
mic web, this happens when gas turns around to form 
sheets (at relatively large densities), and in the spheri- 
cal collapse picture it happens at turnaround (A = 5.55). 
Because our fitting formula in equation (jlOp breaks down 
at that point anyway, our formalism is not meant to de- 
scribe the dense gas most subject to shocks. 

By z ~ 1, shocks become relatively common, and by 
z they dominate the thermal ene rgy of the IGM 
(again see Fig. 11 in iDave et al.l [l999f ) . In the latter 
case, they visibly broaden the scatter around the equa- 
tion of state even at A < 1, refiecting the growing 
prevalence of relatively weak shocks in this low density 
phase. In this epoch, a proper treatment of shocks would 
clearly be necessary to understand the IGM. These weak 
shocks would probably also exist at z > 2, but struc- 
ture formation is much farther behind at that point - 
and flows generally have small enough infall speeds that 
they are unable to shock the gas above the relatively 
high temperatures alread y provided by photoionization 
(IFurlanetto fc Lo"ebll2004f l. In particular, the simulations 
of Dave ct al. ( 19991 ) show that over the redshifts of in- 
terest, gravitational shocks cause negligible scatter in the 
equation of state for A < 10 for the case of hydrogen 
reionization only. This would only hold with greater ef- 
fect once the higher sound speeds due to helium reion- 
ization are taken into account. 

There are also practical difficulties in incorporating 
shocks in our approach. Because the process is stochas- 
tic (in that some, but not all, gas elements of a given 
density are affected), it would require an extra layer of 
sophistication for our model. Those few analytic mod- 
els that do attempt to describe shocks (all of which are 
untested at high redshifts) do not properly account for 
this stochasticity at small d ensities. For examp l e, the 
Press-Schechter approach of IFurlanetto fc Loea (|2004l ) 
associates shocks with gas that has broken off from the 
cosmic expansion and thus sets a density threshold below 
which shocks do not occur. Thus, we will ignore them in 
our calculations, although they may add a small amount 
of scatter (in the form of high temperature gas) to the 
temperature distributions that we show, particularly at 
z = 2; they should be almost negligible at z > 3. 

One additional complication is thermal injection from 
galaxies, such as galactic-scale winds from star forma- 
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tion. These hot shells can both clear out gas and heat 
it. Fortunately, these winds probably only affect regions 
near to galaxies, which are relatively dense (and likely to 
be affected by the nearby ionizing source anyway). We 
therefore neglect these in our calculation as well. 

4. THE TEMPERATURE DISTRIBUTION 

We are now in a position to compute temperature dis- 
tributions of IGM gas elements. For a density (5° , we 
use equation ([5]) to determine the distribution of local 
reionization redshifts and (given zjj and zhc for 
the Universe at large); we then evolve each packet to 
a specified redshift to calculate dp/dT, the temperature 
distribution. 

There is, however, one more ambiguity in our ap- 
proach: each gas element actually has two reionization 
redshifts. To calculate our temperature distributions, we 
must therefore specify some relation between z^ and z^^ 
for a given gas element. One can imagine two simple sce- 
narios. First, the two could be completely uncorrelated. 
That is, a gas element's redshift of hydrogen reioniza- 
tion has no effect on the helium epoch. Alternately, the 
two could be completely correlated: an element that has 
Psp{zion > z^) — 0.1 (or that is ionized quite early in 
the process) will also have Ps {zion > ^hc) ~ ^-^ (also 
early in the process). 

There are strong physical grounds to prefer the latter 
scenario, and we will use it throughout our calculations. 
Intuitively, both redshifts should depend on the nearby 
halo population, and we would expect a region with an 
overabundance of nearby halos at some high redshift to 
have an overabundance of quasars later on. More for- 
mally, we are considering a single mass element, so the 
excursion set trajectory is fixed independent of redshift. 
The reionization barrier is relatively independent of red- 
shift at a fixed ionized fraction (Furlanetto et al. 2006a) , 
so we would expect the trajectory to cross the hydrogen 
and helium barriers at nearly the same ionized fraction. 

Of course, this argument is imperfect because the ex- 
cursion set approach does not properly include stochastic 
fluctuations in the halo population (which are particu- 
larly important for helium reionization). In particular, 
if the number of ionizing sources per discrete bubble is 
less than a few, Poisson fluctuations in that number are 
more important than clustering in determining t he dis- 
tribution of bubbles (see discussion in Furlanctto fc OhI 
I2OO7I) . But it should be a reasonable first guess for our 
analytic model. In reality, the choice has relatively little 
impact on our results, because so long as 2;hc ^ ^^h the 
gas elements approach their thermal asymptotes and the 
precise value of z^ is unimportant (Hui & Haiman 2003) . 

Figure |4] shows the differential and cumulative proba- 
bility distributions for the IGM temperature in pixels at 
the mean density for a variety of redshifts: z — 2, 3, 4, 5, 
and 6 for the solid, long-dashed, short-dashed, dot- 
dashed, and dotted curves, respectively. We have as- 
sumed zh = 8 and Zfje — 3. After helium reionization 
is complete, at z < 3, the distribution has a relatively 
simple form, with a relatively flat distribution between 
T ~ 10^ K and somewhat higher temperatures. This 
is a direct result of the broad distribution seen in Fig- 
ure [TJ most elements were ionized relatively early and 
have cooled to near the thermal asymptote by z — 3. 
Only those elements that were ionized relatively late re- 
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Fig. 4. — Differential and cumulative probability distributions 
(panels a and ft, respectively) of the IGM temperature for a patch 
with (5^ = 0. The solid, long-dashed, short-dashed, dot-dashed, 
and dotted curves are for 2 = 2, 3, 4, 5, and 6, respectively. The 
plot assumes zjj = 8 and zjjo = 3. 



main warm. Both the minimum and maximum of this 
distribution decrease at z < ZHe because of adiabatic 
cooling, until they eventually fill only a small tempera- 
ture range around the thermal asymptote. 

Before helium reionization is complete, however, the 
distribution is more complex. This is largely because 
some elements have undergone helium reionization while 
others have not. The sharp spikes at T ^ 10^ K in panel 
(a) represent the gas that has not yet had its helium 
ionized and so is still approaching the thermal asymp- 
tote from hydrogen reionization. The second component, 
at much higher temperatures and with a similar shape 
to the z = 3 distribution, accounts for elements with 
•^Hc ^ z. Note that this part of the distribution does not 
overlap with the thermal asymptote, especially at the 
higher redshifts. This gas has undergone helium reion- 
ization relatively recently, and cooled since then; hence, 
the minimum of this component decreases with redshift. 
The maximum remains pinned near T ~ TH<= because it 
represents gas that has just been reionized. 

Of course, these distributions also depend on In 
Figure O we show several examples before and after the 
end of helium reionization (at z = 4 and z = 2, respec- 
tively). Again, after zjjc the distributions are relatively 
easy to interpret: regardless of density, the gas has been 
cooling for a substantial interval and so lies near the ther- 
mal asymptote. Note that the densest gas (the solid line 
in Fig. [5?)) is also the warmest at this stage. We have 
already seen that this gas was reionized first, but it is 
also collapsing, and the accompanying adiabatic heating 
has raised its temperature well above that of gas at the 
mean density. 

Again, the distributions are much different before zho : 
a steep rise at low temperatures representing gas with 
z^j, < 4 and a second, smoothly distributed component 
at higher temperatures with z^^ > 4. The relative frac- 
tions in each component are strong functions of density: 
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Fig. 5. — Cumulative probability distributions of the IGM tem- 
perature for patches at 2 = 4 and 2 (panels a and b, respectively). 
In both panels, the curves take 5^ = —5, —2.5, —1, 0, 1, and 2.5, 
from left to right. The plot assumes = 8 and zjjc = 3. In 
the curves in panel a, the sharp rises at T < 10'' K correspond to 
regions that have not yet undergone helium reionization. 



< 20% of the most underdense gas has been ionized 
by 2; = 4, but > 60% of dense gas has been. Again, 
the characteristic temperature of the cool component in- 
creases with density because of adiabatic contraction and 
expansion. The distribution of the fully-ionized, high- 
temperature component is more complex, because the 
adiabatic heating in dense gas competes with the earlier 
ionization redshifts. 

5. THE EQUATION OF STATE 

We have seen above that the competition between adi- 
abatic cooling (or heating) and photoheating implies that 
the IGM temperature depends on the local density, even 
in the absence of inhomogeneous reionization. This is 
usually parameterized by an effective equation of state 



T = To(l + S) 



7-1 



(13) 



where Tq is the temperature at the mean density and 
7 is usually taken to be constant over a broad range of 
density at a given redshift (which turns out to be an 
excellent a pproximation for homogeneous reionization; 
iHui fc Gne din 1997). In that case, the scatter around 
this relation is small, so the equation of state has become 
a primary tool for interpreting the Lya forest. Usually 
To and 7 are taken to be free parameters and fit to the 
observations (although often with priors). 

Inhomogeneous reionization modifies the equation of 
state in two important ways. First, it produces a sub- 
stantial scatter in the temperature, even at a fixed den- 
sity. This implies that a single equation of state will not 
be as accurate a representation of the IGM. Second, the 
slope 7 will change because different densities will (on 
average) be ionized at different times. We will examine 
the latter effect first. 

5.1. Inhomogeneous Reionization and 7 



Figure [6] contrasts the equation of state in the excur- 
sion set model (left panels) with a density-independent 
reionization model (right panels). We assume 0h = 8, 
3, and T^" = 3 x 10'' K throughout. In the top 



panels, we plot the median temperature as a function of 
density.^ In the bottom panels we show the local loga- 
rithmic slope (of the median). 



1 



dlnT 
din (5 ■ 



(14) 



The solid, long-dashed, short-dashed, dot-dashed, and 
dotted curves are for z = 2, 3, 3.5, 4, and 5, respectively. 

First consider density-independent (or stochastic) 
reionization (at right). Here the temperature evolution 
is easy to understand: T is smallest at the highest red- 
shifts, because helium is still only singly-ionized. Once 
helium reionization begins in earnest, the temperature 
increases everywhere by a factor of several. Because 
the post-reionization temperature is independent of lo- 
cal density, this flattens the equation of state, decreasing 
7. In this particular case, the median temperature dis- 
tributions at z — 3 and z = 3.5 are quite similar, but 
that is coincidental. After helium reionization ends, gas 
begins to cool at all densities, but most significantly in 
underdense regions. This translates into an increase in 7. 
At any given redshift, the slope varies only slowly over 
a large range in density. Thus, the assumption of an 
approximately power-law equation of state (as in the ho- 
mogeneous reionization scenario) is probably still valid. 
The main difference here from the usual model is the sub- 
stantial scatter about the median temperature-density 
relation (see §5.2p . 

However, the picture changes completely when the den- 
sity dependence of reionization is included. The most ob- 
vious change is the sharp increase in the median tempera- 
ture as a function of density ed z — 3.5-4. This feature is 
a direct and inevitable result of "inside-out" reionization, 
where ionized bubbles grow around galaxies (or quasars), 
which are highly biased and so include the majority of 
dense gas elements. Specifically, the jump occurs at the 
density at which at least half of the elements have al- 
ready undergone helium reionization. Thus it moves to 
smaller densities as more of the helium is reionized. 

This jump is, of course, accompanied by a discontinu- 
ity in the local logarithmic slope. (This discontinuity is 
a result of choosing the median temperature rather than 
an average; the equation of state defined from the av- 
erage temperature has a much smoother transition but 
similar endpoints at low and high densities.) This point 
also separates two regimes for 7: at smaller densities, the 
curves are fairly steep, because the gas has been cooling 
since hydrogen reionization. Thus adiabatic expansion 
of the voids has steepened the profile. To the right of 
the discontinuity, 7 — 1 « 0.2; this gas has been ion- 
ized relatively recently, with the temperature "reset" to 
a density-independent value, so the slope is smaller. 

Interestingly, density-dependent reionization has 
strong effects even after helium reionization is complete 
(see the z = 2 and z = 3 curves) : 7 — 1 < at relatively 

^ We choose the median because it better matches the methods 
used to mea sure the equation of state directly from Lya forest 
features fe.g.. lSchave et al.l20dOl ). However, the mean may be more 
relevant for some applications, like the power spectrum. 
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Fig. 6. — Equation of state of the IGM. The solid, long-dashed, short-dashed, dot-dashed, and dotted curves are for 2 = 2, 3, 3.5, 4, 
and 5, respectively. For each density, we plot the corresponding median temperature [top panels) and the local logarithmic slope {bottom 
panels). Left: Excursion set model. Right: Assumes that the reionization history is independent of local density (see text). All assume 
2h = 8 and zho = 3, as well as a post-helium reionization temperature of T^'^ = 30, OOOK. 



low densities.^ This inverted equation of state implies 
that voids are hotter than might be expected. This 
region even extends to (5 ~ at z = zhc- This is because 
underdense regions are ionized later in our model and 
have had relatively little time to cool. This trend will 
of course decay after enough time passes, but it persists 
at least until z = 2 in our model, prime territory for 
Lya forest measurements. Moreover it is obvious that 
in this regime it is not appropriate to choose a constant 
value for the logarithmic slope over any substantial 
density interval. Instead, the combination of adiabatic 
cooling, photoheating, and inhomogeneous reionization 
inevitably leads to a complex equation of state. 

Figure [7] shows how the equation of state changes 
with T^p; in the left (right) panels we take T^" = 
4 (2) X 10"* K. Not surprisingly, the initial temperature 
sets the overall scale of the helium reionization break. On 
the other hand, the qualitative features are unchanged, 
and, even in the = 2 x 10'^ K case, we still have 
sharp breaks in the distribution of the median tempera- 
ture, an inverted equation of state at low densities, and 
rapid evolution in the temperature. Artifacts of helium 
reionization remain significant even at z = 2. 

Another uncertainty is that the distribution of de- 
pends on the assumed mass in each gas element, be- 
cause that determines how correlated the (small) mass 
elements are with the (large) ionized bubbles (see Fig. [3]). 
The left panel of Figure[8]shows how this affects our equa- 
tions of state. Although the detailed relations do change, 
all of our qualitative points remain valid. In particular, 
the sharp jumps in the median temperature as a function 
of density still appear, although their precise locations 
change, and 7 — 1 < for small densities, albeit over a 
somewhat smaller range of 5 and z. 

The right panel of Figure [6] shows the evolution in our 

* After zho, the mean and median temperatures are very similar, 
and this negative slope in underdense regions appears in both. 



fast reionization model (with m,„i,i — lOm^, and which 
probably better matches observations). In this case, the 
curves take z = 2, 3, 3.33, 3.65, and 4.35; before zhg, 
these show the same ionized fractions as in Figure[6] (0.7, 
0.5, and 0.2). Again, although the locations of the jumps 
move, the same qualitative features are evident. The 
main difference from our fiducial model is simply that 
the evolution occurs more quickly. 

5.2. Scatter in the Equation of State 

In addition to the change in slope, our model pre- 
dicts a significant increase in the scatter around the me- 
dian equation of state (even excluding the modest scatter 
caused by differenc es in the tidal fields surrounding the 
mass elements; see iHui fc Gnedi3ll997l ). We illustrate 
the scatter at a sequence of redshifts in our fiducial model 
in Figure[9](z = 2, 3, 3.5, and 4, counter-clockwise from 
top left). In each panel, the curves show the equations 
of state for gas with P{< T\5) = 0.1, 0.25, 0.5, 0.75, and 
0.9, from top to bottom. (Thus the solid curves are iden- 
tical to those in Fig. [6l) As may be seen, the scatter 
is substantial and may hamper efforts to determine the 
mean equation of state accurately. 

After helium reionization is complete, the temperature 
scatter in underdense regions is a factor of two or so, 
decreasing toward overdense regions. The scatter is more 
significant during helium reionization. Even if most gas 
elements are on one of the thermal asymptotes, there 
will be significant outliers in the temperature distribution 
(differing by up to an order of magnitude from the mean) : 
in voids, there is a long tail toward high temperatures 
from the rare gas elements with z^^ > z, and in dense 
regions there is a long tail toward low temperatures from 
ga s elements with z^„ < z. Thus, as first pointed out 
bv'H ui fc Haimaiil ([2003) , the scatter in the temperature 
distribution can be used to identify the epoch of helium 
reionization. We have shown here that the density range 
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Fig. 8. — Left: Identical to the left panel of Fig.|6l except we use T = lO"^ K to evaluate the Jeans mass of the IGM. Right: Identical to 
the left panel of Fig. (6] except we use our fast reionization model. The curves are at z = 2, 3, 3.33, 3.65, and 4.35 here; before zjje, these 
show the same ionized fractions as in Figure [6] 



over which this scatter appears wiU evolve from high to 
low densities. This can be used, at least in principle, to 
trace the history of reionization and to understand how 
"inside-out" the process is. 

Of course, the duration over which the scatter per- 
sists will depend on the duration of reionization: for 
example, in the high-mass source model shown in Fig- 
ure [21 the equation of state remains tightly defined until 
z < 3.7. The post-reionization scatter in the high-mass 
source model is smaller, because the range of reionization 
redshifts (and hence cooling times) is smaller. However, 
the differences are not particularly large, and much more 
careful modeling would be required for such quantitative 
constraints to be trusted. In the stochastic model, the 
scatter is also comparable to that in our fiducial model; 



because all gas elements more or less track mean density 
gas, the fractional scatter can be approximately read off 
from Figure [HI Of course, in this case there is no density 
dependence, so there is no jump in the median. 

5.3. Comparison to Previous Work 

To date, the only other theoretical model that allowed 
for inhomogeneous helium reionization is the Monte 
Carlo, semi-analytic model of lGleser et al.l (|2005f ). They 
determined the reionization redshift of volume elements 
probabilistically with the help of a simple biasing pre- 
scription (which, like our model, caused dense pixels to 
be ionized first). They then followed the subsequent evo- 
lution in a similar manner to us, although they allowed 
the local sources to turn on and off and used the log- 
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Fig. 9. — Equation of state of the IGM. In each panel, we show the equation of state taking temperatures from fixed contours of 
10%, 25%, 50%, 75%, and 90% in P(< r|(5° ), from bottom to top. The panels are for z = 2, 3, 3.5, and 4, counter-clockwise from the top 
left. All assume zjj = 8 and zjjo = 3. 



normal model for the density evolution (rather than our 
spherical collapse model). Their fiducial model had a 
slightly faster reionization history than our high-mass 
s ource mo del but only took AT = 10** K. 

iGleser et al.l (120051 ) found a similar overall trend for 
the equation of state (their Fig. 6): it broke into two 
branches, with some high-density pixels that are ionized 
early splitting off to high temperatures and voids staying 
cool until the end of reionization. As in our models, the 
distribution narrows well past reionization. The most 
obvious d ifference is in the a ctual slope of the equation 
of state: IGleser et al.l (j2005[ ) only allowed gas to heat 
by a modest amount during helium reionization, and the 
equation of state never approached isothermality except 
over limited ranges of density. However, Figure [7] shows 
that even in this "weak" heating case, our model gives 
a fairly substantial region of near isothermality (in fact 
even an inverted equation of state). The source of this 
difference is likel y in the prescription for reionization red- 
shifts, for which Icieser et all (j2005( ) used a Monte Carlo 
implementation with a relatively weak bias. Their model 
is thus similar to our density-independent reionization 
scenario, seen in the right panel of Figure [SI where we 
also find no turnover and only modest flattening. 

We conclude that our qualitative conclusions, that he- 
lium reionization leads to a complex and evolving equa- 
tion of state with substantial scatter at low densities, are 
generic to any "inside-out" model of reionization. How- 
ever, the "inversion" in the equation of state is more sen- 
sitive to the detailed prescription for the ionizing sources. 

6. THE TEMPERATURE EVOLUTION DURING HELIUM 
REIONIZATION 

We next turn to the evolution of the IGM temperature 
throughout helium reionization. Of course, we have al- 
ready seen that there is a good deal of scatter in T, and 
it depends on the density of the underlying element. We 
will follow convention and examine the temperature at a 
fixed density (the mean); such regions can be identified 



from their column density in Lya forest spectra (|Schavel 
120011 ) and so constitute a well-defined observational sam- 
ple. 

Figure [TOj shows the temperature as a function of 
redshift in the neighborhood of helium reionization (at 
ZHo = 3). As usual, we take T^'' ^ 3 x 10^ K here. Of 
course, even at a fixed density there is still a range of 
temperatures in our model; we therefore show the me- 
dian temperature at each redshift with the thick solid 
curve, while the thin solid curves show the 10th, 25th, 
75th, and 90th percentiles of the temperature distribu- 
tion. 

The curves split into three regimes. At high redshifts, 
they all lie close together and decrease slowly; this rep- 
resents gas that still has neutral helium. Next there is a 
rapid rise (which, in the case of the thick curve, occurs 
when 50% of the gas has been ionized). The initial spike 
reaches the temperature of the gas that had its helium 
ionized first - because more recent gas elements are now 
hotter. The curves continue to rise, albeit more gently, 
until zhcj because they steadily shift toward gas that 
was ionized later. Finally, at z < zhc, the curves decline 
again, as they simply track a single gas element (ionized 
halfway through the process for the thick curve) while 
it cools toward the thermal asymptote. Note that the 
scatter decreases again in this regime. 

The temperature begins to increase before zhc, with 
the initial spike preceding it by a substantial amount: in 
the case of the median temperature, it occurs at z ^ 4 in 
this particular model. The "reionization redshift" itself 
is best measured by looking at the minimum tempera- 
ture (i.e., the last gas element to be ionized). However, 
in that case the temperature jump is also the least dra- 
matic, because the statistic jumps from gas with neu- 
tral helium to gas that was ionized long before (and has 
cooled significantly) - in this case, AT ~ 6000 K for the 
10th percentile, but AT ~ 10^ K for the median. 

The dotted curve in Figure [TOj shows the median tem- 
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Fig. 10. — Temperature evolution near helium reionization (as- 
sumed to be at 2hc = 3). The thick solid curve shows the me- 
dian temperature for gas elements at the mean density; the thin 
solid curves show the 10th, 25th, 75th, and 90th percentiles of 
the temperature distribution. The dotted curve shows the median 
temperature for a gas element with 5^ = —2.5. 



perature for gas elements with 61 = —2.5 (corresponding 
to a deep void). Because voids are ionized only at the 
end of reionization, the temperature increases closer to 
0Ho, and the jump is somewhat larger (mostly because 
the initial temperature is smaller). 

Figure [Tl] shows how the temperature evolution de- 
pends on the details of the reionization model. In prin- 
ciple, the magnitude of the jump - which increases for 
more rapid reionization - can be used to constrain the 
duration of reionization, but in practice the differences 
may be too small to allow robust inference. In panel (a), 
the solid curves correspond to density-independent reion- 
ization; the dotted curves show the fiducial model from 
Figure \W\ (for the latter, we only show the 10th, 50th, 
and 90th percentiles). The differences are quite modest 
because of course mean-density gas elements have histo- 
ries typical of the universe as a whole, unless reionization 
occurs in a pathological manner. The only obvious dif- 
ference is in the 90th percentile curve at high redshifts. 
With density-dependent reionization, all of the early ion- 
izations occur in dense gas, so these mean density pixels 
get started later. Distinguishing reionization driven by 
density fluctuations from "stochastic" reionization will 
require studying a broader range of densities. 

The bottom panel is similar, except that the solid 
curves use our "fast" reionization model (with mmin = 
lOirii) with the usual density-dependent implementation. 
Comparison to the dotted curves shows that the the tem- 
perature spike occurs at lower redshift (for obvious rea- 
sons - though note that when the 10th percentile is cho- 
sen, the delay almost vanishes) and is also somewhat 
stronger; this is because the median post-reionization gas 
element has undergone less cooling since it was initially 
ionized. The temperature difference is nearly indepen- 
dent of which percentile is chosen. 

Figure [T^ shows the temperature evolution if Tf^° = 
4 X 10** K, again in comparison to the fiducial model. Of 
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Fig. 11. — Temperature evolution near helium reionization (as- 
sumed to be at 2hc = 3). (a) The thick solid curve shows the 
median temperature assuming that the reionization history is in- 
dependent of density; the thin solid curves show the 10th, 25th, 
75th, and 90th percentiles of the temperature distribution, (b) 
Same, except that the reionization history depends on density but 
takes mmin = lOrrii. In both panels, the dotted curves show the 
10th, 50th, and 90th percentile points of the distribution for our 
fiducial model incorporating inhomogencous reionization. 



course, this increases the magnitude of the temperature 
spike without moving its redshift. However, the change 
in temperature is now sensitive to the chosen percentile: 
for example, the peak of the median changes by ~ 6000 K 
while the peak of the 10th percentile changes by only 
~ 2000 K. This is because the 10th percentile curve 
jumps to gas that has already approached the thermal 
asymptote, washing out T^'^. 

Finally, Figure [13] shows how the mean temperature 
(T) evolves over this range (again for gas at the mean 
density). The thick solid curve shows our fiducial model. 
The dotted curves take T/^<= = 4 (2) x lO"* K (top and 
bottom, respectively), the long-dashed curve assumes 
density-independent reionization, and the short-dashed 
curve is for our fast reionization model with mmin = 
lOm^. For comparison, the thin solid curve shows the 
median temperature in the fiducial model. Although 
the mean and the median match almost exactly when 
z < zho, the mean evolves much more smoothly at higher 
redshifts. This is simply because the jump in the median 
occurs when precisely 50% of the gas has undergone he- 
lium reionization, whereas the mean accounts for both 
cold and reionized gas at all times. Interestingly, the 
mean only changes by ~ 9000 K in our fiducial model 
(even though the temperature jump in each gas element 
is > 2 X 10"^ K), and it does so over a fairly large redshift 
interval. 

The differences in (T) between the models are modest; 
the fast reionization model has slightly higher tempera- 
tures after zhc (because less cooling has occurred) and 
slightly smaller temperatures at higher redshifts (because 
less gas has been ionized). The density-independent 
reionization model has a slightly higher temperature at 
high redshifts, because reionization begins a bit earlier 
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Fig. 12. — Temperature evolution near helium reionization 
(assumed to be at Zfjo = 3). The thick solid curve shows the 
median temperature for gas elements at the mean density, with 
= 4 X 10"' K; the thin solid curves show the 10th, 25th, 75th, 
and 90th percentiles of the temperature distribution. The dotted 
curves show the 10th, 50th, and 90th percentile points of the dis- 
tribution for our fiducial model with TH<= = 3 x lO^* K. 
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Fig. 13. — Evolution of the mean temperature (T) near helium 
reionization (assumed to be at ^He = 3). The thick solid curve 
shows our fiducial model. The dotted curve takes Tfo =4x10** K, 
the long-dashed curve assumes density-independent reionization, 
and the short-dashed curves is for our fast reionization model with 
™min = lOrrii. For comparison, the thin solid curves shows the 
median temperature in the fiducial model. 



for mean density pixels. Changing the initial tempera- 
ture obviously changes (T) , though only by ~ 60% of the 
nominal temperature change because of the rapid cooling 
experienced by elements that are reionized early. Inter- 
estingly, if the energy injection is modest, the observable 
change in the mean will be quite small. 

7. COMPARISON TO OBSERVATIONS 

The temperature distribution and IGM equation of 
state are by no means easy to measure, but there have 



been several attempts in recent years. These focus on 
z ~ 2-4, just when our models predict that the distribu- 
tion is most interesting. 

One approach is to measure the temperature directly 
through Lya forest features. Because line widths also de- 
pend on the Hubble expansion, turbulence, and peculiar 
velocities, it is difficult to interpret the velocity width 
attributed to individual absorbers precisely. However, 
one can search for the minimum Doppler parameter for 
systems of a given colum n density and use t hat to mea- 
sure thermal prope rties (jSchave et al.lll999( ). In detail, 
ISchave et al.l ()200Clf ) compared the observed velocity dis- 
tribution with a sample drawn from a simulation that 
did not directly include a model for helium reionization 
but in which the IGM temperature Tq could be varied.^ 
They then measured Tq from the cutoff in the Doppler 
width distribution, using it to gauge the combined ef- 
fects of thermal broadening and Jeans smoothing (the 
latter was estimated using a suite of numerical simula- 
tions). When compared to our model, this process would 
ideally measure the minimum temperature at each red- 
shift (because their simulation did not include the scatter 
from inhomogeneous reionization), although with a finite 
sample it probably actually measures a somewhat higher 
percentile of the distribution. 

Using a s a mple of high-resolution quasar spectra, 
ISchave et al.l ()200Cl( ) found a sharp increase in To, the 
fiducial temperature of mean density gas, at z ~ 3.2: 
from To w 12, 500 K to To ~ 25, 000 K. Interest- 
ingly, the measurements are consistent with a con- 
stant temperature at higher redshifts but a rapidly de- 
creasing temperature at lower redshifts (returning to 
To = 12, 500 K at z = 2); there also appears to be 
some weak evidence that the jump occurs over a finite 
(Az = 0.3) interval. This is certainly qualitatively con- 
sistent wit h helium reionization ending at ZHe = 3.2 
( Schave et al. 2000; iTheuns et al.ll20d2d) . iRicotti et all 
( 2000f) and .McDonald et al.l ()2001[ ) have made similar 
measurements. All are consistent within the errors, al- 
though the latter favored a picture without a sharp jump. 
However, they did not include Jeans smoothing effects 
and also chose lines with well-behaved thermal profiles, 
which may bias the temperature measurement high if 
there is scatter in the equation of state. 

However, there a re a number of su btle difhculties. 
First, the method of lSchave et al.l ()1999D searches for the 
minimum velocity width and attributes it to a combi- 
nation of thermal broadening and Jeans smoothing. We 
would therefore expect that these measurements corre- 
spond to a change in the minimum temperature at any 
given redshift (or at least a low but non-zero percentile) . 
In that case, we expect a smaller temperature jump than 
the median over an extremely short redshift interval. Ac- 
cording to our fiducial model, the observed magnitude of 
the jump (ATo « 12, 500 K) requires Tp - 6, 4, or 
3 X 10^ K if the measurement corresponds to the 10th, 
25th, or 50th percentile, respectively (or slightly smaller 
if reionization is faster). The former two temperatures 
exceed the limit of what might plausibly be achieved even 
with extremely optically thick photoheating. Obviously, 

^ See ITheuns et aP l|2002al ') for a later version that did include 
uniform helium reionization and matched the measurements rea- 
sonably well. 
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we must understand the velocity cutoff in the presence of 
inhomogeneous reionization more precisely, which must 
be done through a careful calibration against simulations. 

A second difficulty is that the observed pre-reionization 
temperature is itself rather large, well above naive pre- 
dictions from models such as ours or from simulations 
(by ~ 4000 K), although the errors are large enough 
that our model cannot be ruled o ut. A similar differ- 
ence occurs in other models (e.g., iTheuns et al.ll200"23 
|Hui_& Haiman 2003'). This discrepancy can be reduced, 
but not eliminated, if hydrogen reionization occurred 
later than z = 8, or it could be due to turbulence, heat- 
ing by a substantially harder spectrum from high-redshift 
quasars, or some other mechanism. We will not try to 
model it in more detail here, deferring a closer look at 
the hydrogen reionization epoch until the future. 

A third difficulty is that the observations also cool 
more quickly than our model at z < zue (albeit again 
at relatively l ow co nfidence), a difference also noted in 
IGleser et al.l (HflOl). This may be a result of evolu- 
tion in the ionizing b ackground, which we have ignored 
(|Theuns et al.ll2002c[) . 

Similar techniques can also be used to measure the 
equation of state of the IGM. So far, the best measure- 
ments lie at z 2-4. Near the mean densitj^it appears 
to be approximately iso thermal at z 3 (jSch avc et aU 
l2000t iRicotti et al.ll2000l ) , with steepening at both higher 
and lowe r redshifts. Interes tingly, the peak in Tq de- 
tected bv lSchave et al.l ()2000[ ) is accompanied by several 
points with best-fit values 7 < 1, although at less than 
la confidence for each one. This is again consistent with 
ZHe = 3.2 according to our model. Because inhomoge- 
neous reionization predicts a complex equation of state, 
with rapid evolution over redshift, it will be interesting 
to study the equation of state over a broader range in 
d ensity. 

IGleser et al.l ()2005() also found that temperature boosts 

larg er than AT = 10* K were necessary to match 

the ISchave et al.l (|2000f ) measurements, although they 
worked with the mean temperature, which we have ar- 
gued is not a good match to the observational methods. 
They also do not reproduce the flattening in the equa- 
tion of state at z 3, which requires a larger temper- 
ature boost or stronger density dependence during the 
reionization process. Overall, it appears that the mea- 
surements require relatively hard sources (or substantial 
IGM filtering) during helium reionization. 

A complementary approach is to use the Lya for- 
est power spectrum of transmitted flux to measure the 
equation of state. Such measurements generally show 
little or no evidence for an evolvi ng equation of state 
or any substantial change in T p (Thcuns et al. l2000l: 
Zaldarri aga et al.ll2001l : IViel et al..,2004 : .McDonald et al l 
2005, 200^, even though they are also most sensitive to 
gas near the mean density. However, the errors are large 
(and in fact the best-fit models often produce unphysical 
results, un less priors are placed on the thermal history; 
IViel et al.l[2004l . 

F igure [T^ helps to res olve this apparent difference with 
the ISchave et al.l (|2000D measurements: the power spec- 
trum is most likely sensitive to the mean temperature, 
rather than the details of the distribution. This evolves 
much more smoothly than the median and undergoes 
a substantially smaller increase near helium reioniza- 



tion (at most, it increases by about half the temper- 
ature jump of an y sing le gas element). For example, 
IZaldarriaga et all (|200lh measure T ~ 2 x 10"^ K at 
z = 3. 9, nearly twice the value estimated bv lSchave et al.l 
(200(fj. However, Figure [13] shows that, in our fiducial 
model, we expect a difference of ~ 5000 K between the 
mean and median (or minimum), which certainly eases 
the discrepancy. The relatively poor redshift resolution 
available from power spectra will further smooth out this 
evolution. More careful modeling is clearly needed, but 
the power spectrum may in fact be consistent with the 
ISchave et al.l (|20M ) measurements . 

Power spectrum measurements probably show greater 
similarity with observations of the evolution in the 
H I effective optical depth re ff(z), which show a 
dip in absorption at z ~ 3.2 (|Bernardi et al.l 120031 : 
iFaucher-Giguere et al.l [20071 ) . presumably due to heat- 
ing from helium reionization (which affects r^g through 
the temperature dependence of the recombination coef- 
ficient, a (X T~°'^). The effective optical depth is also 
more sensitive to the (flux-weighted) mean temperature, 
and appears to require a jump in IGM t emperature at 
heliurn reionization by a factor of ~ 2 (jTheuns et al.l 
l2002aD . even though t^s only changes by ~ 10%. How- 
ever, damping on small scales in the power spectrum, 
which might be expected after an increase in IGM tem- 
perature due to increased Jeans smoothing, is not seen. 
Reconciling these two measurements will require careful 
study with numerical simulations. 

8. DISCUSSION 

We have examined how inhomogeneous reionization af- 
fects the equation of state of the IGM. Because reioniza- 
tion is accompanied by substantial photoheating, it af- 
fects the thermal state long after its completion. We have 
considered models in which the reionization process pro- 
ceeds from large-scale overdensities to large-scale voids 
and in which reionization occurs independently of the 
underlying density. In the former case, dense pockets of 
gas are likely to lie near galaxies (and quasars) and so to 
be ionized relatively early; voids, on the other hand, are 
likely to be the last elements to be ionized. Our semi- 
analytic approach ignores shock-heating in the IGM, as 
well as galactic winds, so it is only accurate at moderate 
or low overdensities (S ^4), where such phenomena are 
relatively rare, but it does describe the bulk of the gas 
contributing to the Lya forest (which sits in photoion- 
ization equilibrium). 

Hydrogen reionization occurred early enough that the 
implications for the thermal state in the observable era 
(z < 4) are relatively modest, and we defer a closer 
look at fluctuations from that process to the future. 
This is primarily because the IGM rapidly approaches 
a (fairly steep) the rmal asymptote, allowing only weak 
const raints on zh (jTheuns et al.l l2002cl : iHui fc Haiman! 
l2003f ). However, He II reionization at z ~ 3-4 has a 
large effect. 

If reionization occurs in a density-independent fashion, 
the temperature always increases with density, but it is 
still nearly isothermal during and near the end of reion- 
ization. If underdense gas is ionized last, the equation 
of state can actually become inverted, with the hottest 
gas elements in the deepest voids. In practice, the "me- 
dian equation of state," defined with reference to the 
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median temperature at each density, is reasonably close 
to isothermal near the end of helium reionization, but 
with a factor of a few spread in the actual temperatures. 
It also develops a break that moves to lower density as 
the process unfolds (although the break is a result of 
taking the median; the "mean" equation of state will 
be smoother). Gas elements denser than the break are 
nearly isothermal (because they have already undergone 
helium reionization) , while those that are less dense have 
temperatures near the thermal asymptote (T < 10^ K). 

As expected, the temperature increases rapidly near 
the end of helium reionization, with the magnitude of 
the jump and its precise timing depending on how it is 
measured. The magnitude of the jump depends on Tj^°, 
but it is significantly smaller than the jump experienced 
by individual gas elements, because some gas elements 
must have been ionized at significantly higher redshift 
and so have cooled toward the new thermal asymptote. 
For many types of observations (for example, those sen- 
sitive to contours in the temperature distribution), the 
observed temperature jump will also occur somewhat be- 
fore reionization is complete. 

Measurements of Lya forest lines do show such a tem- 
perature jump at z w 3.2, accompanied b y a trend to- 
ward a more isothermal equation of state (jRicotti et al.l 
[2000: Schaye et al. 2000). This is qualitatively consistent 
with ZHe = 3.2 in our model, but the magnitude of the 
jump is somewhat higher than one might expect if the 
ISchave et all (|2000( l method measures the minimum IGM 
temperature. If so, we would require an extremely hard 
ionizing background that can produce TP^° ~ 4 x 10^ K. 
It is possible, however, that a more realistic assessment of 
the method with numerical simulations that include in- 
homogeneous reionization would show that it tends to 
measure temperatures closer to the median values, in 
which case a large but not unreasonable initial tempera- 
ture would be required (Tj^° ~ 3 x 10"* K). 

On the other hand, the temperature can also be mea- 
sured from the Lya forest power spectrum, which shows 
no evidence for an evolving equation of st ate an d tends 
to fin d higher temperatures at z ^ 4 (IZald arriae a et al.l 
[Ml lViel et~aII[200l [McDonald et al1l2005Ll200i) . The 
discrepancy probably arises because the power spectrum 
method is sensitive to the mean temperature, which 
evolves more slowly with redshift and has a significantly 
smaller overall jump than statistics like the median. 
Thus inferences about the temperature evolution and 
equation of state will depend on precisely how a mea- 
surement is made (if, for example, the minimum, median, 
or mean temperature at each density is measured, and 
which density range is selected). Careful comparison to 
observations will probably require detailed simulations of 
the reionization process. 

Of course, our model is far from perfect. We have ig- 
nored fiuctuations in the ionizing background (in both 
its amplitude and spectrum), which will create an even 
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larger spread in the gas temperatures. Such fiuctuations 
may also depend on the local gas density and introduce 
biases into the equation of state. For example, voids lie 
far from the ionizing sources, so photons that reach them 
may have been strongly filtered and hence hardened by 
the IGM; on the other hand, after reionization voids ap- 
pear to show systematical ly softer spectra than dense fil- 
aments (jShuU et al.ll2004D . We have also ignored X-ray 
heating before a gas element undergoes helium reioniza- 
tion. In reality, X-rays can travel large distances through 
the IGM and deposit a fraction of their energy as heat. 
This will lift gas off of the thermal asymptote by at least a 
small amount, especially because such photons will tend 
to have higher energies than average, and consequently 
decrease the magnitude of the temperature jump from 
full reionization. It may help to reconcile the overall dis- 
crepancy of our post-hydrogen reionization temperatures 
(at z > 4) with the observations. For example, if a frac- 
tion xhoIII of the helium is fully ionized by photons with 
mean excess energies (above the photoionization edge at 
Eueii) of {E), we would have 

AT^15QQk(^)(4^). (15) 

V 0.1 / V^HcIl/ 

Before helium reionization is complete, this "simmering" 
background is likely to be quite hard (because the soft 
photons would be absorbed near the sources), so the heat 
injection could be substantial. 

There are, at least in principle, a variety of other 
ways to test our model. One could search for density- 
dependent temperature fiuctuations in the Lya forest 
lines, for example through correlations between small- 
scale and large-scale structure in Lya forest lines; the 
small-scale structure becomes a proxy for temperature 
through Jeans smoothing. Such searches have so far been 
unsuccessful (|Zaldarriagal l2002t ITheuns et all I2002bi ). 
but larger samples at a range of redshifts may show a 
signature. The substantial increase in the temperature 
variance during helium reioni zation w ill probably be the 
easiest feature to see (see also lHui fc Haiman 2003), and 
can provide a diagnostic of when helium reionization oc- 
curs and, if the density dependence can be measured, 
how it proceeds. Finally, the equation of state also af- 
fects the power spectrum of the Lya forest, both directly 
(because the optical depth r oc (52-0.7(7-1)^ possibly 
indirectly through spatial flu ctuations in the tempera- 
ture f although .Lai et al.ll2006l have shown that the latter 
are likely to be small). 
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